Singular perturbation analysis of a reduced model for collective motion: 

A renormalization group approach 
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In a system of noisy self-propelled particles with interactions that favor directional alignment, 
collective motion will appear if the density of particles is beyond a critical density. Starting with 
a reduced model for collective motion, we determine how the critical density depends on the form 
of the initial perturbation. Specifically, we employ a renormalization-group improved perturbative 
method to analyze the model equations, and show analytically, up to first order in the perturbation 
parameter, how the critical density is modified by the strength of the initial angular perturbation 
in the system. 
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I. INTRODUCTION 

The interesting phenomena of flocking in animals 
[1] and self-organized patterns in motile cells [^-[1] are 
currently driving the intense theoretical study of collec- 
tive motion among self-propelled particles |9l42]| . Models 
for collective motion usually involve motile particles that 
possess alignment interactions and angular noise. Col- 
lective motion is then observed if the density of particles 
increases beyond a certain threshold. We have previously 
argued that besides density fluctuations, the initial fluc- 
tuations in the heading directions of the particles consti- 
tute another important aspect of the system . Here, 
we determine the critical density as a function of the ini- 
tial perturbation strength by analyzing a reduced model 
for collective motion. Specifically, we assume that the 
angular noise strength, e, is small and employ it as the 
perturbation parameter. We then find that the solution 
obtained by the naive perturbation method is plagued by 
divergences due to the appearances of temporal secular 
terms 23, 24], which we subsequently eliminate up to 
order e by the renormalization group method [25l430l | . 



II. MODEL 

In this work, we consider a minimal model for col- 
lective motion in two dimensions introduced in [20| . In 
this model, every particle is assumed to have constant 
speed u, and that their interactions consist only of a di- 
rectional alignment mechanism, with interaction strength 
£7. Noise, of strength e, is incorporated in the direction 
of travel. Note that the factor e in front of the interaction 
strength is to emphasize that our perturbative treatment 
will be on both the interaction and noise strengths. Since 
7 is still a free parameter, these two effects can be varied 
independently. 



FIG. 1: (a) The form of the initial perturbation with 6 = 0.1 
and rj = 0.5. (b) The form of the function R{x) (c.f. Eq. 
(|22|) '). (c) & (d) The temporal evolutions of a and P, with 
e — 0.5 and 7 = 0.9, respectively. These results are obtained 
by performing numerical integrations of Eqs ([1} with the ini- 
tial square wave perturbation approximated by the function; 
fe[tanh(30(a; + rj)) - tanh(30(a; - j?))]. 
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Starting from the Fokker-Planck equation describing 
this model, we have previously argued that at the onset 
of collective motion, the model equations can be approx- 
imated by a finite set of coupled partial differential equa- 
tions (PDEs) [l^l . Here, we will study the simplest set of 
those coupled PDEs, which corresponds to the following 
two coupled PDEs: 

dta — —2udx(i , dtP — —udxa + e{ja(3 — /3) . (1) 

In the above equations, a corresponds to the local parti- 
cle density, and /3 corresponds to the local vectorial order 
parameter of the system, i.e., a non-zero /3 implies the 
existence of collective motion. Note that the problem is 
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reduced to one dimension as we assume that the initial 
directional preference of the system is along the x-axis. 

The derivation of and the approximations involved in 
the model equations are elaborated in [22^ and we sum- 
marize the essential steps in Appendix for complete- 
ness. We are primarily interested in the regime where 
e <C M and 7 ~ 0{1). The first condition allows us to 
study the reduced model perturbatively with respect to 
e, and the second condition is our main interest because it 
is the region where the transition from disordered motion 
to collective motion occurs, as we will see shortly. 

At the mean field level, i.e., if the spatial variations in 
a and /? are ignored, the threshold for collective motion 
is (c.f. Eqs 0) 0,[23: 
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(2) 



where pc denotes the critical density. We have previously 
argued that such a mean-field picture is incomplete be- 
cause the initial perturbation to the system should play 
a major role as well H^]. For instance, let us assume 
that the initial density is p, i.e., a{x) — p, and 7 < p~^. 
In other words, no collective motion would be expected 
according to the mean-field description. Let us now con- 
sider an initial perturbation in the form of a square wave 
in /3 with magnitude 2b. As shown in the next section, 
such a perturbation would induce a density wave of mag- 
nitude a/26 in a that travels in the positive x direction 
(c.f. Fig. [Tfc) and (d)). Within the traveling density 
wave, a ^ p + \/2b. If 



V2b > 
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(3) 



then the density within the traveling wave is beyond the 
collective motion threshold. Therefore, according to the 
mean-field criterion, we would naively expect that such 
a density wave will be amplified, and thus signals the 
onset of collective motion. In other words, the critical 
density may depend on the strength of the initial per- 
turbation perturbation. In the next two sections, we will 
verify this expectation by analyzing the reduced model 
perturbatively. 



III. NAIVE PERTURBATIVE TREATMENT 

We are concerned with the progression of a pertur- 
bation to an initially disordered system, and so we are 
primarily interested in the initial conditions of the form: 



a{x, t = 0)= A , /3{x, t = 0) = 2 fix) 



(4) 



where A corresponds to the local density and / is an 
arbitrary function that describes the initial perturbation 
in the angular directions of the particles. 

We will now consider the model perturbatively by mak- 
ing the following expansions for a and /3: 



a 



(5) 
(6) 



To the zeroth order, Eqs ([T]) give 

dtao = -2udxl3o , dtPo 



(7) 



and the solutions are the standard traveling wave equa- 
tions: 



ao(x, t) = A + V2 f i^x- \^utj - f {x + V2ut 
j3Q{x, t) = f (^x- V2ut^ + f (x + V2ut^ . 



(8) 



Note that since the units of time and length can be de- 
fined arbitrarily, we will use these degrees of freedom to 
set u — l/y/2 and A = 1 from now on. 
In the first order in e, Eqs ([1]) lead to 

dtai = -V2d^(3i (9) 

dtPi = -^d.ai + 7«o/?o - A • (10) 
v2 

We can eliminate the (3i term above by combining the 
partial derivative of Eq. ([9]) with respect to t and the 
partial derivative of Eq. ([TU|) with respect to x. We 
then obtain: 



(52-a^)ai - y2a,(/3o-7«o/3o) 



(11) 



The above inhomogeneous wave equation is supple- 
mented by the initial conditions; 

ai(x,t = 0) = (12) 
dtai{x,t^Q) = -V29^/3i(a;,t = 0) = (13) 



where the second equality in Eq. p3)) follows from Eq. (|9]) 
and the fact that f3i{x,t = 0) = 0. The general solution 
for ai can thus be written as (e.g., see ch. 8.2 in [3l|): 

ai = 4s / / ' [a,'(/3o ~ -faoM]dx'dt' . (14) 

V2 Jo Jx-(t-t') 

By plugging the expressions for and /3o (c.f. Eqs ([5])) 
into the above integral, we find 

«i = ^[(7-l)(/--/+)+7(./^+./i)] (15) 
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where f± = f{x±t), h± = h{x ±t) and the function h 
is defined by the following ordinary differential equation 
(ODE): 

dyhiy) - fiy) . (16) 

Also, we can deduce an expression for /3i by combining 
Eq. (dni) and Eq. dH): 

/3i = l[(.l-W-+f+)+l{f--fl)] (17) 
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(5- - 9+) 



(18) 



where g± = g(x ±t) and the function g is defined by the 
following ODE: 



9yg{y) = f{y) 



(19) 
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Example 

Consider the case where the initial perturbation is a 
square wave of width 2ri and of magnitude 2&, i.e., 



(3{x, t = 0) = 2b[H{x + 77) - H{x ~ 77)] , 



(20) 



where H{y) is the Heaviside function such that H{y) — 1 
for y > 0, and zero otherwise. Note that in order for the 
perturbative treatment to be vahd, the magnitude b has 
to satisfy the fohowing conditions: eb,dP <C u = l/-\/2 
(c.f. Eq. (H))). For this square wave perturbation, we have 
from Eqs and ([T! 



h± = b^Rix ± t) , g± = bR{x ± t) . (21) 
where (c.f. Fig. [Tfb)) 



, y <-ri 

R{y) = { y ^ -'V<y <v 
'2'n , y>v- 



(22) 



The temporal evolutions of a and /3 as a result of this 
initial perturbation are shown in Fig. [Ijc) and (d). 



IV. RENORMALIZATION GROUP 
TREATMENT 

The appearances of t in front of the terms in squared 
brackets in Eqs ([T5|) and ([T7|) signify that the pertur- 
bative solution only makes sense for t <C e~^. In other 
words, the problem is a singular perturbation problem 
and the terms in squared brackets are called secular terms 
[23 |. To eliminate such secular terms, we will now follow 
the renormalization group (RG) method introduced in 

MM- 

In the RG method, we first introduce an arbitrary time 
T and split i as i — r + r. We then rewrite f{xzLt) as 
B±{r)f± where /± = f{x ± {t - C{t)), such that 



1 = 



l + ^5±(r)e" 



71=1 



B±{r) 



= C{t)^Y.''^^^>" ■ 



(23) 
(24) 



n=l 



The series {6^} and {c„} are the multiplicative and ad- 
ditive renormalization constants respectively, and they 
will be chosen order by order in e to eliminate the terms 
containing r. 

Applying the above expansions to a and ignoring the 
non-secular term for the time being, we have for up to 
order e, 

a = l + y/2[{l + h^e)B-f--{l + hte)B+f+] (25) 
"7-1 



^/2 



■{B.f.-B+f+) 



(26) 



+7(B-,/-+B+/+ 
7-1 



V2 



B-f- - B+U) + 7(B-./- + 5+,/;) 



(27) 



(28) 



Our strategy now is to eliminate the term shown in Eq. 
(1^51) by picking b^ appropriately. We will again focus on 
the square wave as our initial perturbation, i.e.. 



/(y) = H{y + ry) - H{y - r?) 



(29) 



The virtue of this choice is the property that = f±- 
It is due to this critical property that we can perform the 
elimination simply by defining the following: 



7-1 
V2 



B± + iBl 



With b^ defined as above, a is simplified to 

a = 1 + \/2[B_/_ - 
"7-1 



+l{Blf_ + Blfl 



The same expansion renders /3 as follows: 



+6(t-r) 



{B^f^+B+f^ 



.J^{Blf_-Blfl) 



(30) 

(31) 
(32) 

(33) 
(34) 
(35) 



The benefit of the above exercise is that as t is arbi- 
trary, the divergence of the original secular terms can 
be eliminated by picking r to be close to t. Now, the 
next step is to also eliminate r in the above equations. 
As a and /3 should be independent of how r is chosen, 
their derivatives with respect to r are zero. Therefore, 
to first order in e, we arrive at the following differential 
equations governing B± and C: 



f- 



f- 



dB^ 



dB^ 
~d7 



dB+ 
~d7 



dB+ 
~d7 



dC 
97 



e 

, 7-1 



^{Blf^+BlU) 
{B^f^-B+f^ 



V2 



V2 
7-1 



+ 



^{Blf^-BlU) 
{B^f^+B+,U 



V2 



. 



(36) 



(37) 



(38) 



Eq. ([55]) indicates that we can simply set C(r) to be zero. 
Furthermore, Eqs ([55]) and (1571) can be decoupled and we 
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arrive at the following two ODEs: 

7-1 





e 


dr 




dB+ 


e 


dr 


" 71 



7^! 



V2 



B4 



The corresponding solutions are 

7-1 



B+{r) 
B-ir) 



7-1 

5_e-E(7-l)r/2 _ 



(39) 
(40) 

(41) 
(42) 



where b± are obtained from the initial condition of the 
problem. 

By picking t — t, we can finally eliminate the original 
secular terms, and the expressions for a and /3 become 



1 



/3 = B-/- + B+l 



2(i?_/_ - B+U) + e^{h^ - h+) + 0{e^) 



-(.9_-5+) + 0(e^) 



Note that in the above solutions, we have reinstalled the 
non-secular terms h± and g±, which are defined in Eqs 
((2T|l . We show that h± and g± are not modified under 
the renormalization procedures in Appendix IB] 

We now consider the temporal evolutions of the system 
under the designated initial perturbation. By assump- 
tion, we have at i = 0, 



a(x, t = 0) = 1 , I3{x, t^0) = 26/(0) 
Hence, B±{Q) = b and so 



7-1 



b- = 



7-1 



+ \/27 



(43) 



(44) 



Let us focus only on the density wave in a traveling to 
the right. The magnitude of the density wave is 



B-{t) 



6(7-1) 



((7 - 1) -t- \/267)e-<^(T-i)*/2 - 7267 



(45) 



If 7 > 1 , then _B_ grows with time since the first term 
in the denominator approaches -^267 as t grows [13) ■ In 
other words, the original perturbation is amplified, which 
is a signature for the onset of collective motion. This 
scenario corresponds exactly to our mean-field prediction 
(c.f. Eq. ©). The more interesting case is that even for 
7 < 1, B- still grows if 



6 > 



1-7 
V27 



(46) 



This is the condition discussed previously and based 
purely intuition ground (c.f. Eq. ([3])). In other words, 
the critical density in the presence of the perturbation 
considered is shifted downwards and becomes 



Pc=l 



V2b . 



(47) 



V. CONCLUSION 

In this paper, we have demonstrated how the critical 
density for collective motion can depend on the strength 
of the initial angular perturbation. Our conclusion is 
based on a renormalization-group improved perturbative 
treatment of a reduced model for collective motion. At 
the simulation level, our results suggest that simulations 
with different implementations of the initial perturba- 
tions may lead to different values of the critical density. 



Appendix A: Reduced model 

For completeness, we summarize here the essential 
steps in deriving the equations describing the reduced 
model studied here [22| . We consider a minimal model 
for collective motion in two dimensions, where every par- 
ticle is assumed to have constant speed, u, and that their 
interactions consist only of a directional alignment mech- 
anism. Noise, of strength e, is incorporated in the direc- 
tion of travel. Specifically, let there be N particles in a 
volume of V, their equations of motion are: 



dr, 
dt 

dt 



2u , , 
-v(e.) 

TT 

9t/(R,e) 



(Al) 
(A2) 



where I < i < N , H = (n, . . . , tat), 6 ee (6*1, . . . , On), 
w{9) = (cos 9, sin 0), and the noise is assumed to be Gaus- 
sian characterized by the following moments: 

(77,(t)) = , {m{t)ri,{t')) - 5,,5{t - t') . (A3) 

Moreover, the alignment interaction is assumed to be of 
very short range and can thus be approximated by a delta 
function: 

UiR, 8) = y - r,) coa{e, - 6^) , (A4) 

l<j 

where £7 corresponds to the coupling strength. If we 
denote the probability distribution of the density of par- 
ticles in the state (R, 6) at time t by /(t,R, 6), then 
the Fokker-Planck equation corresponding to the system 
is 1331: 



9/ 
dt 



d 



Focusing now on the single-particle density function, p, 
where 



P(ri,6'i 



(TV!) / drz ■ ■ ■ drjvd^z • • • dgjv/(R, 9) 
(A^-1)! 



(A6) 
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Eq. (jA5l) becomes 



as follow: 



2u 

TV 



cose ^ + sm^ ' 



dx 



dy 



TV de 



(A7) 



(jV!)/dr3---drjvdg3---dgjv/(R,e) 
(TV -2)! 



dt 



where 

p(2)(ri,^i,r2,(?2) = 

. (A8) 

We now employ the molecular chaos assumption to 
close the above hierarchical equation by approximating 
p^'^\r,9,r,9') by p{r,6)p{r,6'). Within this approxima- 
tion, the Fourier transform of Eq. (|A7[) with respect to 
the angular variable, 9, leads to 

dtPn{r) = -en'^pnir) - u dx{pn+i{r) + pn-i{r)) 
+idy{pn-i{r) ~ /5„+i(r)) 



/5-i(r)/5„+i(r) - pi(r)p„_i(r) (A9) 



where p(r, 9) = J2'^=~oo Pnir)e '"^ and /5„(r) is complex. 
In [i^], we have argued that at the onset of collective 
motion, only the lower modes are important. Therefore, 
we truncate the above set of infinitely many equations by 
ignoring all p„ such that n > 1. Furthermore, if we pick 
the initial angular perturbation to be directed towards 
the positive x direction, then we need only consider the 
real part of p+ i , as the imaginary part of p±i will never 
be excited [22|. With these simplifications, the original 
single particle density function is approximated as 



p{r,9) ~ a{r) + 2;3(r)cos( 



(AlO) 



where a(r) corresponds to the local density, and /3{r) 
corresponds to the local vectorial order parameter since 



(v(r)) - (( 



cos f , sm 



(/3(r),0) 



(All) 



From Eq. (|A9I) , the reduced model for the onset of collec- 
tive motion is therefore described by the following equa- 
tions: 

dta = -2ud^l3 , dtp = -ud^a + e(7a/3 - /?) . (A12) 



a 



V2 



h + Oie^) (Bl) 

P = B^f^ + B+}+ + + 0{e^) (B2) 

where h,g are functions of x,t to be determined by Eqs 
(P). We note that d'^B{t) = 0(e"), and so to order e, 
we have from Eqs ^ the following differential equations 
governing h,g: 

= V2{E-f- - E+f+) + ^dth + -^d.g (B3) 



= - E+f+ + -dtg + -d^h (B4) 

-e[(7 - l){B-f- + B+f+) + V2(B2 /_ - BlU)] 



where E± = e^^dtB±. Similar to the derivation of Eq. 
([M]) . the above equations can be combined to give the 
following second order differential equation for h: 



(B5) 



-2e 



(7-l)(i?-/:+B+/;) (B6) 



2{BUL - Bif^) 



= 



(B7) 
(B8) 



where the last identity follows from the definition of E± 
(c.f. Eqs ^ and ^). Given that h{x,t ^ 0) = and 
dth{x,t = 0) = -2(S_/_ - E+f+) (c.f. Eq. (|B3l)), we 
have 

h{x,t) = 2 r^\E+{Q)-E_mf{s)ds (B9) 

Jx-t 

rx+t 



-26^7 / f{s)ds 

Jx-t 

l{h- - h+) , 



(BIO) 
(Bll) 



where / is as defined in Eq. (|29p . and h± are as defined 
in Eqs (j2ip . Moreover, given the expression for /i, we can 
now employ Eq. (|B3|) to deduce that 



g{x,t) = (1 - 7)(g- - g+) 
where g± are also defined in Eqs (HI]). 



(B12) 



Appendix B: Non-secular terms 
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